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^^ ■ An algorithm is proposed for constructing a group (ensemble) pulsar time based 

bJQ. on the application of optimal Wiener filters. This algorithm makes it possible to 

■^ \ separate the contributions of variations of the atomic time scale and of the pulsar 

^^ • rotation to bary centric residual deviations of the pulse arrival times. The method 

is applied to observations of the pulsars PSR B1855+09 and PSR B1937+21, and is 

used to obtain corrections to UTC relative to the group pulsar time PTgns- Direct 

Mh' comparison of the terrestial time TT(BIPM06) and the group pulsar time PT, 

O ' 

Vh ■ shows that they disagree by no more than 0.4 it 0.17 /is. Based on the fractional 

C^ , instability of the time difference TT(BIPM06) - PTens> a new limit for the energy 

— 1 \ density of the gravitational- wave background is established at the level Q^gh^ ~ 10~^. 
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The discovery of pulsars in 1967 [1] and of millisecond pulsars in 1982 [2], as well as 
subsequent observations, have clearly shown that the rotational stability of pulsars is such 
O \ that pulsars can be used as astronomical clocks. Currently, the accuracy with which pulse 

arrival times (PATs) can be measured is at the level of several microseconds for most pulsars, 
. and reaches the submicrosecond level for some pulsars. If we compare this accuracy to the 

H ' interval covered by observations, we obtain the relative accuracy 10~^^, which is essentially 

comparable to the fractional instability of the best atomic frequency standards. 

A number of studies have been concerned with the stability of pulsar rotation and its 
application to time systems. Writing about the principles for establishing Terrestrial Time 
(TT), Guinot [3] concludes that Terrestrial Time as realized by the Bureau International de 
Folds et Mesures [TT(BIPMXX), where XX denotes the year] is the most expedient time 
system to use for pulsar chronometry. The principles of pulsar time and the definition of the 
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pulsar second are given in [4]. It is shown in [5] that it is not possible to use the rotation of 
pulsars to improve the definition of the time unit. The studies [6-9] introduce a definition 
for a time based on the motion of a pulsar in a binary system - Binary Pulsar Time, or 
BPT - together with theoretical expressions for the variations in the rotational and orbital 
parameters as functions of the duration of the observing interval. The main conclusion of 
these studies is that BPT is less stable than the ordinary pulsar time PT over short intervals, 
but the fractional instability of BPT can reach 10"^'^ over long intervals (10^ -^ 10^ yr). Petit 
and Tavella [10] present an algorithm for determining a group pulsar time based on weighted 
averaging, which is compared to an algorithm based on optimal filtration; they also consider 
some ideas connected with the stability of BPT. Foster and Backer [11] present a polynomial 
approach to describing clock and ephemerides variations and the infiuence of gravitational 
waves passing through the solar system on pulsar chronometry. 

Here, we propose a method for forming a group pulsar time (PT) based on the applica- 
tion of optimal Wiener filters [12, 13]. Section 2 discusses the pulsar-chronometry (timing) 
algorithm from the point of view of time scales. Section 3 describes the principles for con- 
structing optimal filters and the formation of a group pulsar time. Section 4 contains the 
results of computer simulations of detecting a signal against the background of correlated 
noise using optimal filtration. Section 5 discusses the results of applying Wiener filters to 
observations of the pulsars PSR B1855+09 and PSR B1937+21. 

1. PULSAR CHRONOMETRY ALGORITHM 

An observer located on the Earth, which is rotating about its axis and revolving around 
the Sun, receives a signal from a pulsar using a radio telescope over an accumulation time 
that is sufficient to obtain a specified signal-to-noise ratio. The pulse arrival time (PAT) 
is measured in a local time scale via crosscorrelation with a standard profile for the pulsar 
pulse. The resulting PATs tn can be translated into UTC, TAI, and TT using the relations 
[14] 

UTC = Tjv + At, TAI = UTC + A;, TT = TAI + 32.184 s, (1) 

where Ar is the local time correction. A; is a whole number of seconds introduced to take 
into account the variation in the length of the day, and the quantity 32.184 s is added to 
remove a historical jump between atomic and ephemerides time. Since the duration of the 



second in TT depends on the position and velocity of the Earth in its orbit, an additional 
transformation from TT to Barycentric Time TB is required, using the algorithm described 
in [15]. 

The observations that have been transformed into the uniform Barycentric Time system 
TB must then be reduced to the barycenter of the solar system, in order to exclude variations 
in the PATs due to the motion of the observer about the barycenter [14] 

T = t-to + Ania, 5, /i«, /x^, tt) + Aorb - DM/f + A,ei, (2) 

where to is the initial epoch; t the topocentric PAT in TB; T the PAT at the barycenter of 
the solar system; A^(q;, 6, fia, f^s, tt) the Remer correction along the Earth's orbit; a, S, Ha, f^s 
and vr the right ascension, declination, proper motion in right ascension and declination, and 
parallax of the pulsar, respectively; Aorb the Remer correction along the orbit of the pulsar, 
when it is located in a multiple system; DM / f'^ the delay for the signal propagation in the 
interstellar and interplanetary media at frequency / taking into account the Doppler shift; 
and Arci the relativistic correction for the delay of the signal propagation in the gravitational 
field of the solar system. 

The PATs at the solar-system barycenter are then used to calculate the rotational phase 
(number of rotations) of the pulsar: 

N{T) = No + yT+hT^ + e{T), (3) 

where N^ is the initial phase at time T = 0, u and z> are the rotational frequency of the pulsar 
and its derivataive at T = 0, and e{T) represents variations in the rotational phase (timing 
noise). The procedure for determining the parameters includes refining a, 6, fia, f^s and vr, via 
a least-squares fit minimizing a weighted sum of the square differences between N{T) and 
the nearest integer. The residual deviations in the rotational phase are usually expressed in 
units of the time 6t = 5N/v. Here, we consider variations in the intrinsic rotation of the 
pulsar and variations due to irregularity of the reference time scale Atciock{T). 

When intercomparing different realizations of atomic time [3], flicker noise in the frequency 
dominate on intervals of several months, while random walk noise in the frequency dominates 
on intervals of several years. Thus, clock variations have a power spectrum of the form l/cj" 
in the frequency domain, and are expressed in the time domain by the polynomial 

Atdock(r) = Co + cT + -cT^ + -dT^ + .... (4) 
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It is clear that the appearance of the quantity Atciock in Eq. (3) leads to a redetermination 
of the rotational parameters of the pulsar: 

N{r) = K+{1 + c)fT + \{fc + (1 + cff)T\ (5) 

where T is ideal barycentric time and / and / are the rotational frequency of the pulsar and 
its derivative undistorted by the clock parameters. For this reason, we must use the PAT 
values expressed in TT(BIPMXX), as best determined at the current time [5]. 

2. OPTIMAL FILTRATION 

Let us consider n uniform measurements of a random quantity (the residual PAT devia- 
tions) r = (ri, r2, . . . , r„). The quantity r is the sum of two uncorrelated quantities, r = s+e, 
where s is a random signal to be estimated and associated with the clock contribution, and 
e is the contribution of the rotational phase of the pulsar, which we treat here like noise. 
We are interested in estimating the signal s against the background of the additive noise e 
using the Wiener filtration method. 

Wiener filtration consists of estimating the signal s given the measurements r and covari- 
ances (7) [16]. We must reconstruct the random signal with insufficient a priori information, 
since the covariance matrix of the signal is not known a priori, and is estimated from the 
observational data themselves via a crosscorrelation of all the data, assuming that the vari- 
ations in the clicks (the estimated signal) and in the pulsar rotational phase (additive noise) 
are uncorrelated. 

The optimal Wiener estimate of the signal s is given by [16]: 

s = Q,,Q;;r = Q,,Q;> = Q,,(Q,, + Q,e)-'r, (6) 

where the covariance matrices Q^r, Qsd and Q^^ are formed as arrays of the corresponding 
covariance functions. The covariance functions for r, s and e are written 

cov(r,r) = {ri,rj) = {si,Sj) + {ei,ej) , 

cov(s, s) = (Sj, So) , 

^ ' ^ '' (^,j = l,2,...,n). (7) 

cov(s,r) = {si^Tj) = {ri,Sj) = {si,Sj) , 

COv(£,e) = {Ei^Ej) . 



The angular brackets, (), denote ensemble averaging. We assume that the processes s and e 
are stationary in a weak sense; i.e., only the first and second moments are stationary, since 
fitting a quadratic polynomial to the rotational phase excludes the nonstationary part of 
the random process [17]. Thus, the covariance functions depend only on the time difference 

Since the pulsar-chronometry algorithm assumes that the PATs are defined relative to a 
reference time, distinguishing the covariances {si,Sj) and {Si^Ej) requires observations of at 
least two pulsars using a single time system. In this case, combining the pulsar PATs and 
assuming that the cross-covariance (^^j, ^Sj) = i^Si^'^ej) = 0, we can estimate 

(s„ sj) = ((V, + \, V, + V,) - (V, - V„ V,- - V,)) /4 (8) 

or 

(s„s,) = (V„V,). (9) 

If M pulsars are observed to construct a group pulsar time, we will have M{M — l)/2 
estimates of the covariance matrices (sj, Sj) = i^ri, Vj) , (A;, / = 1, 2, . . . , M) 

The matrix Q",.^ in (6) is a whitening filter. The matrix Q^^ forms a signal from the 
whitened data. 

The averaged signal (group pulsar time) is written 

M(M-l) 

2 ~^ ^ 



M{M 



^ Y. „Q..-5^^«;U";-'r, (10) 



m=l j=l 

where *w is the relative weight of the ith pulsar, *w = k/Uj^, 0"j is the rms deviation of the 
whitened data *Q~^^ -^r, and the constant k serves to ensure that ^^ % = 1. The first factor 
in (10) is the average cross-covariance function, and the second factor is the weighted sum 
of the whitened data. 

The following algorithm was used to calculate the auto- and cross-covariances. The ob- 
servational data ^Ti were subjected to a rapid Fourier transform, 

n 

^x(a;) = ^Y. '^thte''^'^ (A; = 1, 2, . . . , M), (11) 

where the weights /i^, which are bell-like functions, were used to reduce leakage through the 
sidelobes [18]. They can be calculated to very good accuracy using the formula [18] 



J„ {i,\V(n - l)\/[l-(2^-lP] 



"' - ^'"^ ^(.H-(»-l)) -' <^'' 
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where Cq are scaling factors used to ensure that ^ /i^ = 1, Jq, is a modified zeroth-order 
Bessel function of the first kind, and the parameter W serves to regularize the sidelobe level 
in the spectrum, and is usually taken to be in the range W = 1 ^ A. We used the value 
W = 1. 

The power spectrum {k = I) and cross-spectrum {k ^ I) were calculated using the formula 

''Xico) = ^\'xicoyx*ico)\, (13) 

where (■)* denotes complex conjugation. 

Finally, the auto- {k = I) and cross- {k ^ I) covariance were calculated using the formula 



U)=l 



^X{uj)e-'''\ {k,l = 1,2,..., M). (14) 



3. COMPUTER MODELING 

We carried out numerical simulations to estimate the quality of the signal reconstruction 
applying Wiener filtration and weighted averaging [10]. In constrast to [13], where a harmonic 
signal was taken as the initial signal, we estimated a time series of the difference UTC- 
TT(BIPM06) in the interval MJD = 46399-48949, from which we subtracted the quadratic 
polynomial from the least-squares fit. We added white noise with zero mean and various 
dispersions. For example, if we used 50 pulsars for the modeling, the dispersion of the 
white noise Hq was o"^ = 1 ^ 50. The weights were taken to be inversely proportional to the 
dispersion of the obtained noise series. The quality of the reconstructed signal was calculated 
as the rms deviation of the difference between the reconstructed and input signals. 

Figures la and lb show the quality of the signal reconstruction using the weighted aver- 
aging (dashed curve) and Wiener filtration (solid curve) for white noise, as functions of the 
number of pulsars used and the length of the data series. As the length of the data series 
and the number of pulsars increase, the quality of the Wiener-filtration signal reconstruc- 
tion becomes increasingly higher compared to the weighted-average reconstruction as the 
rotational phase variations increase. 

4. RESULTS 

The method described above was applied to observational data (barycentric residual de- 
viations of the PATs) for the pulsars PSR B1855+09 and PSR B1937+21 [19]. Since these 
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Figure 1. The rms difference between the reconstructed and initial signals as a function of (a) the 
number of pulsars M and (b) the length of the data series n, using weighted averaging (dashed) and 

Wiener filtration (solid). 

data are not uniform, we applied a spline interpolation to regularize the data, with the aim 
of simplifying the subsequent reduction using matrix- algebra methods. We chose a step of 10 
days, which preserved the original quantity of data. This transformation of the data distorts 
the high-frequency component, but leaves the low- frequency component that is of interest 
to us unchanged [20]. 

We took the part of the data that was common to both pulsars (251 points for each), in 
the interval MJD = 46450-48950, for this reduction. Since, generally speaking, the series of 
residual deviations have different means and slopes after their lengths are shortened, they 
were reduced by fitting a quadratic polynomial (Fig. 2). 

The data for PSR B1855+09 and PSR B1937+21 were obtained in UTC based on the 
information from [19]. Consequently, the signals that were distinguished from the observa- 
tional data for PSR B1855+09 and PSR B1937+21 were the corrections UTC-PT1855 and 
UTC-PT1937 . The combined signal (group time) UTC-PTens is shown in Fig. 3 by the 
thin curve. It shows behavior similar to that for UTC-TT(BIPM06), with the correlation 
coefficient p = 0.75 ± 0.04. A direct comparision of terrestrial time TT(BIPM06) and the 
group Pulsar Time PTgns shows that they disagree by no more than 0.4 ± 0.17 fis. 

We calculated the fractional instability cr^ of the time difference TT-PTens, which was 
az = (0.5 ± 2) ■ 10^^^ on a time interval of seven years. 
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Figure 2. Barycentric residual deviations ^'^r (in /is) for PSR B1855+09 (points) and PSR 
B1937+21 (curve) after fitting a quadratic polynomial. 
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Figure 3. UTC-TT(BIPM06) (in /is, bold curve) and UTC-PTons (thin curve). 



5. DISCUSSION 



The fractional instability of time systems is characterized by the so-called Allan variance, 
which is numerically equal to the second-order finite difference of the clock phase variations. 
A chronometric analysis of observational data contains a determination of the rotational 
parameters of the pulsar, including at least the first frequency derivative; this corresponds 
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to fitting a quadratic polynomial in time to the rotational phase of the pulsar, which, in 
turn, is equivalent to excluding the second-order derivative from the PAT data. For this 
reason, the use of the classical Allan variance is not expedient. Instead, another quantity 
was proposed as a means of characterizing the fractional instability of the rotation of pulsars 
- (Jz [21]. A detailed algorithm for calculating (Jz is given in [22]. 

Figure 4 plots the fractional instability of the intrinsic rotation of the pulsars PSR 
B1855+09 and PSR B1937+21 taking into account the contribution of the reference time 
and the variations TT-PTcns- Theoretical curves for the behavior of a^ when noise due to 
the presence of a stochastic gravitational-wave background with relative energy densities 
VtgK^ = 10"^ and 10^^" [19] are also shown in the lowerright corner of the figure. The az 
curve crosses the line Qgh"^ = 10^^°; however, the upper limit for the gravitational-wave 
background should be an order of magnitude higher, taking into account the uncertainty in 
this quantity. Thus, as a conservative estimate of this upper limit, we adopt Qgh!^ < 10^^. 

The fractional instability of TT(BIPM06) relative to PTcns is at the level < 10~^^ for a 
seven-year interval, and is an order of magnitude better than the fractional instability of the 
rotation of PSR B1855+09 and PSR B1937+21, taking into account the contribution of the 
reference time system. As the numerical simulations show, when Wiener filtration is used, 
the accuracy with which TT-PTcns is estimated grows with the number of pulsars used. We 
estimate the accuracy of the optimal filtration method to be 170 ns. This accuracy was 
calculated as the rms deviation between the obtained and smoothed signals. The smoothing 
was carried out in the frequency domain using a low-frequency Kaiser filter with a bandwidth 
of /max/32, where /max = 2/At and At = 10 days is the sample interval. Our choice of this 
bandwidth was based on the desire to obtain smoothness similar to the UTC-TT(BIPM06) 
curve. Thus, the uncertainty in our estimation of PTcns can, in principle, reach several tens 
of nanoseconds if the most stable millisecond pulsars are used. 

The proposed method does not distinguish quadratic trends in the reference time and in 
the rotational phase of the pulsar, because the pulsar itself displays secular deceleration of the 
rotational frequency. In our opinion, this does not hinder the use of pulsars as independent 
clocks, since longer-period variations in the reference clocks will be revealed as data are 
accumulated over increasingly longer time intervals. 

The low relative accuracy of P noted in [5] likewise does not pose problems, since the 
rotational phase of the pulsar is not predicted. However, if it is required to predict the 
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Figure 4. Fractional instability cr^, based on the chronometry data for PSR B1855+09 (solid) and 

PSR B1937+21 (dashed) and corresponding to TT-PTons (dot-dashed) as a function of the 
observing interval t (in years). The theoretical curves of cr^ for flgh^ ~ 10^^ and 10^^° are shown in 

the lower-right corner of the figure. 



behavior of some reference atomic time, such as UTC, relative to the group pulsar time, this 
can be done based on variations of UTC-PTcns using standard methods for the prognosis 
of time series, provided that this prognosis is obtained relative to short-period fluctuations, 
without linear and quadratic trends. In this approach, the low relative accuracy of the 
derivative of the rotational period of the pulsar will not play a role, since the absolute phase 
is not predicted. 
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6. CONCLUSION 

We have presented an algorithm for forming a group pulsar time based on optimal Wiener 
filtration. The algorithm makes it possible to distinguish the contributions to barycentric 
residual deviations of PATs due to irregularities in the intrinsic rotation of the pulsar and 
variations in the reference clocks. Both irregularity in the pulsar rotation and variations in 
the reference time scale are obtained relative to an ideal time system. Realization of the 
algorithm requires observations of at least two pulsars relative to a common reference time. 

The proposed approach has better accuracy than the use of weighted averaging to form 
the group time scale, since it uses additional information about the signal via its correlation 
function or power spectrum. 

The presence of a pulsar time that is independent of terrestrial conditions makes it pos- 
sible to carry out independent checks of terrestrial time scales, which immediately provides 
possibilities for revealing the presence of variations that are common to all terrestrial stan- 
dards that would otherwise be undetectable. 
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